/****

This routine picks a massive star from the Kroupa (2001) IMF, which is 
   dN / dlog m \propto 
    m^0.7 (0.01 < m < 0.08)
    m^-0.8 (0.08 < m < 0.5)
    m^-1.7 (0.5 < m < 1)
    m^-1.3 (1 < m < 120)
 
We truncate the IMF to m > 0.08 since there are no models for that regime.

In case you want to extend 0.01 < m < 0.08, you can adopt these weights
  //#define WGT1 0.4967
  //#define WGT2 0.4360
  //#define WGT3 0.0426
  //#define WGT4 0.0247


Written by MF&RD, Apr 2010

*******/

//these are the weights for m > 0.08

#define WGT1 0.       
#define WGT2 0.866
#define WGT3 0.0846
#define WGT4 0.049

double kroupa_imf(gsl_rng *rand) {
  
  double x, mmin, mmax, gamma;

  //First decide which power-law interval to draw from. Note we use
  //gamma 1 less than the exponent in Kroupa (2001), because we want
  //to draw from dN/dm = (dN / dlog m) / m
 
  
  x = gsl_rng_uniform (rand);
  
  if (x <= WGT1) {
    mmin = 0.01;
    mmax = 0.08;
    gamma = -0.3;
  } else if (x <= WGT1+WGT2) {
    mmin = 0.08;
    mmax = 0.5;
    gamma = -1.8;
  } else if (x <= WGT1+WGT2+WGT3) {
    mmin = 0.5;
    mmax = 1.0;
    gamma = -2.7;
  } else {
    mmin = 1.0;
    mmax = 120.0;
    gamma = -2.3;
  }

  // Now draw from a power-law in m and return 
  return(ranpowerlaw(mmin, mmax, gamma, rand));
}
